rm(list=ls(all=TRUE))

path <- "" #here you need to adjust the path macro to the place on your computer
setwd(path)

library(haven)
library(dplyr)
library(ggplot2)
library(Hmisc)
library(lme4)
library(texreg)
library(effects)
library(estimatr)
library(sjPlot)
library(lmtest)
library(plm)

LevNum <- function(x){
  c(round(quantile(x, na.rm=T)[2], digits=5),
    round(mean(x, na.rm=T), digits = 5),
    round(quantile(x, na.rm=T)[4], digits=5))
}

dat <- read_dta(paste0(path, "\\party_governments.dta"))

dat$iso2code <- as_factor(dat$iso2code)
dat$country <- as_factor(dat$country)

# adding some new variables

dat$mainstream[dat$family2==2 | dat$family2==4 | dat$family2==6 | dat$family2==8] <- 1
dat$mainstream[is.na(dat$mainstream) & !is.na(dat$family2)] <- 0
dat$mainstream <- as.factor(dat$mainstream)
table(dat$mainstream, useNA ="always")

dat$past_pm[dat$pm_counter>=1] <- 1
dat$past_pm[dat$pm_counter<1 | is.na(dat$pm_counter)] <- 0
dat$past_pm <- as.factor(dat$past_pm)
table(dat$past_pm, useNA ="always")

dat$opp_party <- as.factor(dat$opp_party)

# low and high contexts

cab_cntr <-  dat %>% 
  group_by(iso2code, cabinet_name) %>%
  summarise_at(vars(all_evts, all_parties_evts), funs(mean(., na.rm=TRUE)))

cab_cntr <- within(cab_cntr, prt_share <- all_parties_evts*100/all_evts)
table(cab_cntr$prt_share)
cab_cntr$cabinet_name[cab_cntr$prt_share>100]
quantile(cab_cntr$prt_share, na.rm=TRUE)[3]

diff.cntry <- c("BE" ,"DK" ,"FI" ,"FR" ,"GR" ,"IE" ,"IT" ,"LU" ,"MT" ,"NL" ,"NO" ,"PT" ,"RO" ,"SI" ,"UK")
nodiff.cntry <- c("AT" ,"BG" ,"CY" ,"CZ" ,"EE" ,"DE" ,"HU" ,"IS" ,"LV" ,"LT" ,"PL" ,"SK" ,"ES" ,"SE" ,"CH")

cab_cntr$prt_share_cat[as.character(cab_cntr$iso2code) %in% diff.cntry] <- 0
cab_cntr$prt_share_cat[as.character(cab_cntr$iso2code) %in% nodiff.cntry] <- 1

cab_cntr$prt_share_cat <- relevel(as.factor(cab_cntr$prt_share_cat), "0")

dat_upper_lev_cat <- subset(cab_cntr, select=c("iso2code",
                                               "cabinet_name",
                                               "prt_share",
                                               "prt_share_cat"))
dat <- merge(dat, dat_upper_lev_cat)
rm(dat_upper_lev_cat, cab_cntr)

table(dat$prt_share_cat)

# Party level regression

##################################
###### Appendix D table 2 ########
##################################

m.1 <- lm(parties_evts ~
              left_right*prt_share_cat +
              mainstream*prt_share_cat + 
              past_pm*prt_share_cat + 
              opp_party*prt_share_cat + 
              vote_bef*prt_share_cat +
              vote_ch_perc_bef*prt_share_cat +
              length_gov + country,
              data=dat)

m.2 <- lm(parties_evts ~ 
            state_market * prt_share_cat + 
            liberty_authority * prt_share_cat  +
            mainstream * prt_share_cat + 
            past_pm * prt_share_cat + 
            opp_party * prt_share_cat + 
            vote_bef * prt_share_cat + 
            vote_ch_perc_bef * prt_share_cat  +
            length_gov + 
            country,
            data=dat)

m.3 <- lm(parties_evts ~ 
            state_market * prt_share_cat + 
            liberty_authority * prt_share_cat  +
            mainstream * prt_share_cat + 
            past_pm * prt_share_cat + 
            opp_party * prt_share_cat + 
            vote_bef * prt_share_cat + 
            vote_ch_perc_bef * prt_share_cat  +
            mass_org_scale * prt_share_cat +
            length_gov + 
            country,
            data=dat)

m.1.cf <- coeftest(m.1, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m.2.cf <- coeftest(m.2, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m.3.cf <- coeftest(m.3, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))

htmlreg(list(m.1, m.2, m.3),
          override.se = list(m.1.cf[,2], m.2.cf[,2], m.3.cf[,2]), 
          override.p=list(m.1.cf[,4], m.2.cf[,4], m.3.cf[,4]), 
          omit.coef="country",
          caption = "",
          single.row=TRUE,
          custom.note="*** p < 0.001, ** p < 0.01, * p < 0.05<br/>All models include country fixed effects.<br />Standard errors clustered by countries.",
          custom.coef.map = list("mainstream1" = "Mainstream",
                                 "past_pm1" = "Past PM party",
                                 "opp_party1" = "Opp. party",
                                 "vote_bef" = "Vote share",
                                 "vote_ch_perc_bef" = "Ch. in vote share",
                                 "left_right" = "Left-right",
                                 "state_market" = "State-market",
                                 "liberty_authority" = "Liberty-authority",
                                 "mass_org_scale" = "Mass party org.",
                                 "length_gov" = "Length of Gov.",
                                 "prt_share_cat1" = "Differentiation",
                                 "prt_share_cat1:mainstream1" = "Mainstream*Diff",
                                 "prt_share_cat1:past_pm1" = "Past PM party*Diff",
                                 "prt_share_cat1:opp_party1" = "Opp. party*Diff",
                                 "prt_share_cat1:vote_bef" = "Vote share*Diff",
                                 "prt_share_cat1:vote_ch_perc_bef" = "Ch. vote share*Diff",
                                 "left_right:prt_share_cat1" = "Left-right*Diff",
                                 "state_market:prt_share_cat1" = "State-market*Diff",
                                 "prt_share_cat1:liberty_authority" = "Liberty-authority*Diff",
                                 "prt_share_cat1:mass_org_scale" = "Mass party org.*Diff"),
        file=paste0(path, "\\prt_reg_inter_NA.html"))

##################################
###### Appendix D table 1 ########
##################################

dat_diff <- dat[dat$prt_share_cat==0,]
dat_nodiff <- dat[dat$prt_share_cat==1,]

m.1d <- lm(parties_evts ~
            left_right +
            mainstream + 
            past_pm + 
            opp_party + 
            vote_bef +
            vote_ch_perc_bef +
            length_gov + country,
          data=dat_diff)

m.2d <- lm(parties_evts ~ 
            state_market + 
            liberty_authority  +
            mainstream + 
            past_pm + 
            opp_party + 
            vote_bef + 
            vote_ch_perc_bef  +
            length_gov + 
            country,
          data=dat_diff)

m.3d <- lm(parties_evts ~ 
            state_market + 
            liberty_authority  +
            mainstream + 
            past_pm + 
            opp_party + 
            vote_bef + 
            vote_ch_perc_bef  +
            mass_org_scale +
            length_gov + 
            country,
          data=dat_diff)

m.1nd <- lm(parties_evts ~
            left_right +
            mainstream + 
            past_pm + 
            opp_party + 
            vote_bef +
            vote_ch_perc_bef +
            length_gov + country,
          data=dat_nodiff)

m.2nd <- lm(parties_evts ~ 
            state_market + 
            liberty_authority  +
            mainstream + 
            past_pm + 
            opp_party + 
            vote_bef + 
            vote_ch_perc_bef  +
            length_gov + 
            country,
          data=dat_nodiff)

m.3nd <- lm(parties_evts ~ 
            state_market + 
            liberty_authority  +
            mainstream + 
            past_pm + 
            opp_party + 
            vote_bef + 
            vote_ch_perc_bef  +
            mass_org_scale +
            length_gov + 
            country,
          data=dat_nodiff)

m.1d.cf <- coeftest(m.1d, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m.2d.cf <- coeftest(m.2d, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m.3d.cf <- coeftest(m.3d, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m.1nd.cf <- coeftest(m.1nd, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m.2nd.cf <- coeftest(m.2nd, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m.3nd.cf <- coeftest(m.3nd, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))


htmlreg(list(m.1d, m.1nd, m.2d, m.2nd, m.3d, m.3nd),
        override.se = list(m.1d.cf[,2], m.1nd.cf[,2], m.2d.cf[,2], m.2nd.cf[,2], m.3d.cf[,2], m.3nd.cf[,2]), 
        override.p=list(m.1d.cf[,4], m.1nd.cf[,4], m.2d.cf[,4], m.2nd.cf[,4], m.3d.cf[,4], m.3nd.cf[,4]), 
        omit.coef="country",
        custom.model.names = c("M1 (High Diff.)", "M1 (Low Diff.)", "M2 (High Diff.)", "M2 (Low Diff.)", "M3 (High Diff.)", "M3 (Low Diff.)"),
        caption = "",
        single.row=TRUE,
        custom.note="*** p < 0.001, ** p < 0.01, * p < 0.05<br/>All models include country fixed effects.<br />Standard errors clustered by countries.",
        custom.coef.map = list("mainstream1" = "Mainstream",
                               "past_pm1" = "Past PM party",
                               "opp_party1" = "Opp. party",
                               "vote_bef" = "Vote share",
                               "vote_ch_perc_bef" = "Ch. in vote share",
                               "left_right" = "Left-right",
                               "state_market" = "State-market",
                               "liberty_authority" = "Liberty-authority",
                               "mass_org_scale" = "Mass party org.",
                               "length_gov" = "Length of Gov."),
        file=paste0(path, "\\prt_reg_split_NA.html"))

